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ABSTRACT 


There  is  an  observed  variation  in  body  wave  magnitude  of  approxi¬ 
mately  —  0.2  units  for  explosions  in  different  areas  of  Yucca  Flat, 
Nevada.  The  variation  appears  to  correlate  with  the  Cenozolc  basin 
structure  at  Yucca  Flat.  The  basin  has  been  modeled  geophysically  by 
Ferguson  (unpublished  manuscript) .  In  this  study  it  is  hypothesized 
that  the  variation  in  m^  is  caused  by  scattering  or  resonance  effects 
within  the  local  geologic  structure.  This  conjecture  has  been  in  - 
vestigated  by  computation  of  synthetic  teleseismic  P-wave  amplitude 
responses  for  the  Yucca  Flat  geophysical  model.  The  technique  due  to 
Aki  and  Lamer  (1970)  was  used.Good  quantitative  agreement  with 
observations  was  found. 
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INTRODUCTION 

Body  wave  magnitude  estimates  for  a  number  of  nuclear  explosions 
at  Yucca  Flat,  Nevada  have  been  analyzed  for  spatial  variation 
(Alewine,  personal  communication) .  Each  explosion  was  corrected  to 
a  common  datum  by  the  use  of  well  established  scaling  laws  for  yield 
and  depth  of  burial  (Mueller  and  Murphy,  1971).  These  data,  when 
plotted  on  the  map  of  Yucca  Flat  in  figure  22,  show  a  progressive 
decline  of  about  O.S  magnitude  units  from  the  middle  of  the  valley 
toward  the  east  over  a  distance  of  four  kilometers.  This  variation 
complicates  the  analysis  and  calibration  of  seismic  yield  estimation 
and  source  discrimination  techniques. 

The  magnitude  variation  is  correlative  with  the  structure 
contours  of  the  top  of  the  Paleozoic  units  as  presented  in  the  FIB 
model  of  Yucca  Flat.  The  systematic  nature  of  the  variation  suggests 
a  deterministic  explanation  of  the  phenomena  in  terms  of  scattering 
of  the  outgoing  body  waves  by  the  local  structure. 

This  problem  is  not  soluhle  in  any  simple  form  for  realistic 
geophysical  models  of  Yucca  Flat.  A  numerical  solution  might  be 
computed  by  any  of  several  techniques  such  as  finite  difference, 
finite  element  or  some  form  of  ray  tracing.  After  a  review  of  the 
literature  on  numerical  wave  propagation  techniques,  the  above 
alternatives  were  rejected  in  favor  of  the  Akl  and  Larner  (1970) 
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Fig.  22.  Contours  of  body  wave  magnitude  variation  (m^)  super¬ 
imposed  on  a  contour  map  of  the  Tertiary-Paleozoic  interface  for 
Yucca  Flat,  Nevada.  The  m^  contours  were  modified  from  Alewlne 
(personal  communication) . 
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collocation  mat hod.  This  technique  proalsed  to  have  high  accuracy 
and  reasonable  computation  speed.  It  is  not  limited  to  plain  strain 
theory  like  the  finite  difference  calculation  or  high  frequency  like 
the  asyntotic  ray  theory;  also  the  solution  is  not  restricted  to  a 
limited  set  of  rays  or  nodes. 


DERIVATION  OF  THE  METHOD 


In  this  section  a  technique  introduced  into  seismology  by 
Akl  end  Lamer  (1970),  Lamer  (1970),  Bouchon  (1973)  and  Bouchon  and 
Aki  (1977),  will  be  developed  in  a  general  form.  The  method  assumes 
a  solution  in  each  homogeneous  layer  or  region  in  terms  of  a  series 
of  plane  waves.  The  boundary  conditions  between  the  regions  are 
matched  at  a  finite  number  of  points  on  the  boundaries.  In  numerical 
analysis  this  type  of  calculation  is  called  an  orthogonal  basis 
boundary  value  collocation  scheme.  The  solution  is  expanded  in  an 
orthogonal  basis  of  known  solutions  to  the  partial  differential 
equation  and  forced  to  satisfy  the  boundary  conditions  at  a  finite 
number  of  points.  The  interpolatory  characteristics  of  the  solution 
are  used  to  provide  an  approximate  satisfaction  of  the  boundary 
conditions  at  unsampled  boundary  points.  Due  to  certain  requirements 
of  the  reciprocal  method  that  will  be  used  to  apply  the  explosion 
source,  it  is  more  convenient  to  develop  the  method  in  terms  of 
potentials.  From  Lame's  theorem  (Aki  and  Richards,  1980,  p  68),  the 
equations  of  motion  for  a  homogeneous  elastic  medium  will  yield 
solutions  in  terms  of  a  scalar  and  a  vector  potential. 


la)  U 
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each  k  and  only  three  scalar  potentials  are  required  (Lamer,  1970, 
p  25-44).  The  new  propagation  coordinates  are: 


30  X  -  U1.V)X. 

3»)  IXx'l  • !  Kx  •»  Tjyj  , 

3e)  K  s  X  COS  O*)  , 

COS  C&) 
sin  (12) 

O 


sin  (si) 
cos  (Ifl) 
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The  displacements  are  written  In  terms  of  the  required  potentials  as: 
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/c)  <f)  •  Va  j  , 

Id)  <j>  *  /$*  V2  f  . 

These  equations  include  separate  wave  equations  for  compressional , 
lc,  and  shear.  Id,  waves.  In  general  the  potentials  are  ambiguously 
determined  and  hence  are  useful  only  when  certain  symmetries  are 
present  in  the  formulation. 

Plane  wave  solutions  with  harmonic  time  dependence, 

i)  <j>  -  A  exp  (-  i  (cM  +  kx  ■*  *)y  -  ^)) } 

will  be  assumed.  In  the  case  of  a  wave  incident  on  a  horizontal 
Interface  all  scattered  wave  vectors  will  be  contained  in  the  same 
vertical  plane  as  the  incident  wave.  An  alternative  statement  is 
that  k  and  T)  are  constant  for  all  waves  satisfying  the  boundary 
conditions.  Thus  the  coordinates  may  be  rotated  to  eliminate  one 
wavenumber  component  and  the  potentials  reduced  to  three  scalar 
potentials,  with  horizontally  polarized  shear  (SH)  motion  uncoupled 
from  compresslonal  (P)  and  vertically  polarized  shear  (SV)  motion 
(Aki  and  Richards,  1980,  p  215).  In  the  case  of  an  interface  with 
topography  variable  in  the  x-direction  only  and  incident  waves  making 
the  angle  SI  with  that  direction,  the  y  wavenumber  component,  f)  ,  will 
remain  constant,  but  the  x  wavenumber  component,  k  ,  will  couple  over 
the  entire  wavenumber  spectral  range.  The  constancy  of  y)  permits  the 
specification  of  a  vertical  plane  for  each  value  of  k  which  will  con¬ 
tain  the  scattered  waves.  The  P-SV  and  SH  motions  will  uncouple  for 
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In  each  homogeneous  region  the  general  solution  will  be  expended  In 
an  Infinite  number  of  plane  wave  contributions: 

5c)  ^  Cx.y.t.k)  *  [  A ,  (k)  exp  (-'»>0  ♦  Az  Cu)  «p  0^*^  «p{*fa+9y))> 

Sb)  *  Cx,y,i»  a  [a,  (u)  «*P  A  (W)« p  (>>*)]  «p(-i(kxt?yi), 

Sc)  X  Cx,y,i,k)  =  |/\3  (k)  «pC-«7»)  *■  A„  (k)cxp  (hri)]  exp(-t(kx^y)). 

the  term,  exp (-lot),  is  suppressed  In  all  the  following  equations. 

The  boundary  conditions  to  be  satisfied  are  continuity  of 
displacement  and  traction  normal  to  the  boundary  at  each  Internal 
Interface.  The  normal  traction  is  set  to  zero  on  the  free  surface. 
The  general  solutions  for  the  displacement  and  stress  components  are: 

U)  U;  (.x.y.it)  *  Ajt  txp  (r'kx)  dK) 

<*b)  O'-  Cx.y,*)  *  j  A*  exp(-.Kx)  dk, 

' .  j  >k  ,m  =  1,1,3,  3  1  )•••>£»  , 

The  normal  tractions  at  an  interface  are  obtained  from  the  inner 
product  of  the  unit  normal  vector  and  the  stress  tensor  or  the 
operator, 

7)  Fii  s  "Fiji,  (I)  >  *  i  j  *  '>*>■*>  -*s  /,2j. 

The  unit  normal  is  a  function  of  the  interface  topography,  |  (x). 

The  operators  e,  f  and  F  are  exponential  functions  of  y,  z,  k  and  ^  , 
as  well  as  the  elastic  constants  and  density.  They  are  wavenumber 
domain  spaclal  derivative  operators  as  required  by  4.  The  internal 
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boundary  conditions  require  that  these  quantities  be  evaluated  on 
each  side  of  the  Interface.  In  matrix  form  this  is  expressed; 
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exp  (-ikx)  dk  . 


The  integration  over  k  must  be  performed  numerically,  so  the  Fourier 
transform  is  replaced  by  a  Fourier  series  of  finite  length  (rectangle 
rule  integration  or  discrete  Fourier  transform) .  This  results  in  a 
periodic  repetition  of  the  interface  topography  with  fundamental 
wavelength,  L.  The  collocation  technique  evaluates  the  displacements 
and  tractions  at  a  finite  number  of  points  in  the  interval  [o,l]. 

The  space  and  spectral  discretizations  are: 

«W)  Xn  =■  Ax  r>  ,  n  =  O,  />•••>  fO  -/  , 

Km  3  m 2, 


At  this  point  the  notation  required  becomes  rather  complex.  In 
addition  to  the  tensor  indices  we  now  have  indices  associated  with 
the  series  summation  and  collocation  points.  The  differing  layers 
of  the  model  must  also  be  distinguished.  Figure  23  is  introduced  to 
clarify  some  of  the  additional  notation.  The  superscripts,  i,  will 
serve  to  identify  Interfaces  starting  with  the  free  surface,  i-1. 

The  superscripts,  1,  refer  to  the  layers,  where  layer  1  is  just  below 
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Fig.  23.  An  Introduction  to  the  indlcial  notation  used  to 
identify  the  homogeneous  regions  and  Interfaces  used  In  the  Akl- 
Larner  technique  derivation. 


12 


the  lch  interface.  There  are  P  layers  and  Interfaces.  For  the 
discrete  space  specified  by  9,  equation  8  nay  be  written, 

io)  U*  *  (*  A' . 

In  10,  U*  is  the  6*N  vector  of  displacements  and  tractions  evaluated 
at  the  discrete  Interface  points,  Gi  is  the  6  •  M  x  6  •  N  matrix 
which  transforms  the  discrete  1c,  6  *  M  vector  of  potential 
coefficients.  A1,  in  layer  1.  In  the  simplest  application  we  will 
set  M-N  for  a  square  linear  system  of  equations.  A  similar  equation 
exists  in  layer  i  +  1  at  the  same  interface,  which  would  be  the 
(i  +  l)th  interface.  If  the  U*-  and  vectors  evaluated  at 
\  (x)  are  set  equal  the  result  may  be  solved  for  A1  in  terms  of 


,0  A'  ■  [§']" 

The  product  matrix  is  the  propagator  matrix  for  this  problem.  The 
product  of  each  such  propagator  at  all  interfaces, 

I  -  &' «')  TT  (V)]'1  O')  , 

-  -  jsq 

results  in  a  relationship  between  the  displacement  and  stress  at  the 
free  surface  and  the  wave  field  in  the  half-space,  layer  P,  given  by 

•i)  u!  *  J  A f 

An  incident  body  wave  is  specified  in  the  half  space,  then  the  free 

surface  boundary  conditions  of  vanishing  stress  are  used  to  solve  for 

P 

the  reflected  waves  in  the  halfspace  thus  completely  specifying  A  . 


The  A*  vectors  la  each  other  layer  are  found  by  application  of  the 
propagator  matrices  in  equation  11. 

Lamer  <1970)  introduced  a  variation  of  the  collocation  method 
which  is  of  considerable  importance.  The  partitions  of  0  and  the 
columns  of  the  partitions  of  G  are  Fourier  transformed  over  N  space 
samples.  The  boundary  conditions  are  then  specified  In  the  wave-* 
number  domain.  The  resulting  Fourier  series  is  then  truncated  to  MSN 
points  and  the  equations  solved  as  before,  except  for  an  inverse 
transform  to  convert  the  displacements  and  tractions  back  to  the 
space  domain.  For  M*N  this  procedure  is  identical  to  the  space 
domain  collocation,  otherwise  it  Is  a  least  squares  solution  of  the 
collocation  equation. 

The  wavenumber  spectrum  for  this  problem  possesses  poles 
corresponding  to  surface  wave  modes.  Integration  along  the  real 
wavenumber  axis  will  result  in  numerical  instability.  The  inclusion 
of  viscoelastic  attenuation  will  remove  the  poles  from  the  real 
k  -  axis  and  permit  integration  along  contours  parallel  to  that  axis. 
The  attenuation  may  be  specified  as  an  average  temporal  Q  for  the 
model  with  a  complex  frequency  and  wavenumber  or  as  an  average 
spatial  Q  with  complex  velocity  and  wavenumber.  In  order  that  the 
proper  phase  relationships  exist  across  the  incident  wavefront  the 
imaginary  part  of  the  horizontal  wavenumber,  in  all  layers,  must 
equal  the  imaginary  part  of  the  horizontal  wavenumber  of  the  incident 
wave  (Aki  and  Lamer,  1970).  The  imaginary  part  of  the  frequency  is 


determined  by. 


•4)  X™  (p)  =  \lrr  If |  /(2Q)  , 

for  a  temporal  Q. 

The  numerical  accuracy  of  Che  solution  technique  ia  very  high.' 
If  the  interface  topography  is  adequately  sampled,  so  that  no 
aliasing  occurs  (i.e.,  the  Fourier  expansions  are  convergent),  then 
the  solutions  are  computed  with  accuracy  equivalent  to  the  rounding 
error  of  the  computer.  The  simultaneous  equations  In  11  and  13  are 
usually  well  conditioned,  so  that  little  loss  of  significance 
occurs  in  their  solution.  A  different  form  of  error  results  from 
the  so  called  Rayleigh  ansatz  (Aki  and  Lamer,  1970).  For  interface 
slopes  greater  than  about  60°,  the  solution  should  contain  multiply 
reflected  waves  which  are  not  provided  for  in  the  specification.  In 
practice  this  error  has  not  been  very  important. 


SOURCE-RECEIVER  RECIPROCITY 


The  Aki-Larner  formulation  may  be  used  to  calculate  the 
responae  of  an  Irregularly  layered  model  to  an  incident  plane  body 
wave.  The  theory  must  be  extended  to  include  point  explo8lve 
sources.  Bouchon  and  Aid.  (1977)  demonstrate  how  a  discrete  wave¬ 
number  representation  can  be  used  to  calculate  the  near  field 
response  of  complex  sources  in  irregularly  layered  structures. 

Because  of  our  interest  in  teleselsmic  measurements  and  the  far  field 
response,  a  more  efficient  procedure  results  from  the  use  of  the 
seismic  reciprocity  theorem  (Bouchon,  1976)  and  the  plane  wave 
response.  The  reciprocity  theorem,  after  White  0965),  states: 

If,  in  a  bounded,  inhomogeneous,  anisotropic  elastic  medium, 
a  transient  force  f(t)  applied  in  some  particular  direction 
at  some  point  P  creates  at  a  second  point  Q  a  transient 
displacement  whose  component  in  some  direction  is  u(t), 
then  the  application  of  the  same  force  f(t)  at  point  Q  in 
direction  will  cause  a  displacement  at  point  P  whose 
component  in  the  direction  is  u(t). 

The  response  of  a  point  dilatatlonal  source,  located  at  P  in 

layer  j,  at  a  point  Q,  in  the  half  space  layer  n,  can  be  derived  from 

the  consideration  of  a  force  P  at  Q.  The  force,  F,  is  directed  along 

the  ray  connecting  P  and  Q  and  results  in  a  vector  displacement  U(P) 

at  P.  By  the  reciprocity  theorem  the  same  force  at  P,  in  the 

displacement  direction,  will  produce  a  displacement  U(Q)  in  the  ray 

direction  (P-wave  motion) .  The  dilatatlonal  source  is  modeled  by 

three  orthogonal  force  dipoles.  To  obtain  the  displacement  at  Q,  in 


16 

ray  direction,  we  compute  the  dilatation  at  P  and  modulate  by  M,  the 
desired  source  moment  spectrum  (von  Seggern  and  Bland ford,  1972). 

«  sc®  ■  f- (  ^  +  + 

Ue  substitute  the  potential  at  P  due  to  a  plane  wave  incident  along 
the  desired  ray  Into  15  to  obtain, 

lb)  U  (9)  *  "ZyT  (A,"*  e*P  exp («!>**))  exp(-'«0<x™>y>). 


Equation  16  provides  the  teleselsmlc  P-wave  motion  due  to  the 
explosive  source.  The  efficiency  of  this  method  derives  from  the 
ability  to  calculate  the  response  from  any  source  location  for  a 
given  ray  with  only  one  model  plane  wave  response  calculation.  In 
contrast,  the  discrete  wavenumber  source  representation  would  require 
a  re-evaluation  of  the  propagator  matrices  for  each  source  location. 
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A  THEORETICAL  RESPONSE  PROFILE 

% 

The  procedure  outlined  in  the  previous  sections  has  been  coded 
in  FORTRAN  IV.  The  program  has  been  tested  on  a  number  of  simple 
problems  Including  flat  layered  and  gently  perturbed  structures. 
Working  versions  exist  for  both  the  Control  Data  Corporation  6600 
and  the  Digital  Equipment  Corporation  VAX  11-780. 

The  scattering  hypothesis  can  be  tested  on  a  profile  taken 
from  the  Yucca  Flat  model.  The  first  calculation  has  been  done  on 
an  east-west  profile  at  a  latitude  of  37°  02'  N.  The  basin  has 
very  little  variation  parallel  to  strike  and  good  structural  control 
is  provided  by  geophysical  surveys  in  that  area.  In  figure  24,  the 
observed  m^  variation  south  of  37°  05'  is  plotted  along  with  a 
profile  taken  from  Goforth,  et.  al.  (1979).  The  calculated,  one 
Hertz  P-vave  amplitude  variation  for  sources  distributed  along  the 
profile  in  figure  24  should  agree  with  this  profile  if  the  hypothesis 
is  correct. 

The  geologic  structure  of  the  profile  can  be  approximated  by 
four  layers.  The  Paleozoic  limestone  can  be  modeled  as  a  half  space 
which  is  faulted  at  a  low  angle  (~50°)  on  the  west  by  the  Carpetbag 
fault  system  so  that  an  asymmetrical  basin  is  formed.  The  basin  is 
filled  with  Tertiary  volcanics  and  Quaternary  alluvium.  The  water 
table  at  a  depth  of  about  500  meters  divides  the  volcanics  into 
saturated  and  dry  units.  The  alluvium  is  entirely  free  of  water. 
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Fig.  24.  A  profile  across  Yucca  Flat,  Nevada  at  a  latitude  of 
37°  02'  N  showing  the  observed  m.  variation.  The  geological  struc¬ 
ture  was  taken  from  Goforth,  et.al  (1979). 
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Average  velocities  and  densities  for  these  units  were  derived  from 
well  logs  as  contained  in  Table  III.  The  profile  is  shown  at  the 
bottom  of  figure  25. 

The  response  due  to  sources  at  the  four  locations  shown  in 
figure  25,  for  P-wave  angles  of  emergence  between  -40  and  40  degrees 
(west  and  east)  and  propagation  vectors  within  the  model  plane,  were 
computed.  An  average  Q  of  50  was  assumed  at  a  frequency  of  one 
Hertz  and  a  source  depth  of  0.8  kilometer  was  utilized.  This 
calculation  resulted  in  a  polar  radiation  pattern  of  amplitude  vs. 
emergence  angle  for  each  shot  point  as  shown  at  the  top  of  figure  25. 

In  order  to  model  the  averaging  process  inherent  in  the 
magnitude  calculation  the  median  amplitude  of  each  shot  will  be 
compared.  The  median  is  plotted  as  a  semi-circular  arc  on  each 
radiation  pattern.  The  ratio  of  median  amplitudes  between  shots  on 
the  west  (near  the  Carpetbag  fault)  and  those  on  the  east  is  2.4. 

This  corresponds  to  an  m^  variation  of  0.4  magnitude  units,  in  good 
agreement  with  the  profile  in  figure  24. 

The  results  of  figure  25  are  limited  in  several  important  ways. 
Ho  out  of  plane  rays  are  considered  and  only  a  single  narrow  frequency 
band  was  studied.  Future  calculations  will  include  the  out  of  plane 
rays  and  a  time  domain  synthesis  over  a  broader  (one  or  two  octave) 
frequency  band.  The  radiation  patterns  display  a  rather  complicated 
variation  which  may  be  observable  in  teleseismlc  observations.  This 
phenomena  can  only  be  studied  by  examination  of  a  number  of  specific 
shots  and  receivers,  with  spectral  and  time  domain  modeling  of  the 
observed  P-waves. 


TABLE  III 


VELOCITY  STRUCTURE  OF  YUCCA  FLAT,  NEVADA 


Stratagraphlc 

Unit 

P-Wave 

Velocity 

S-Wave 

Velocity 

Poisson's 

Ratio 

Density 

QAL 

1.34 

0.64 

0.35 

1.80 

TV  (Dry) 

2.14 

1.14 

0.30 

1.80 

TV  (Wet) 

3.00 

1.60 

0.30 

1.80 

PZ 

4.57 

2.64 

0.25 

2.50 

Units  are  km/sec  and  gm/cc. 


wmm 
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Fig.  25.  Polar  radiation  diagrams  of  teleseismic  P-wave  amplitude 
as  a  function  of  emergence  angle  for  shots  distributed  across  Yucca 
Flat,  Nevada. 
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CONCLUSIONS 


The  results  of  this  investigation  indicate  that  it  is  possible 
to  model  the  effects  of  near-source  scattering  on  teleseisalc  signals 
from  Yucca  Flat,  Nevada.  This  is  made  possible  by  the  existence  of  a 
reasonably  accurate  geophysical  model  of  the  geologic  structure 
developed  from  gravity,  borehole  and  seismic  exploration.  The  Aki- 
Larner  technique  is  a  computationally  effective  means  of  performing 
response  calculations  for  quite  complicated  models,  such  as  the  one 
exhibited  here.  Future  efforts  will  concentrate  on  time  domain  wave¬ 
form  synthesis  and  detailed  source  receiver  pair  models.  It  is  also 
possible  to  model  other  profiles  across  the  valley.  Similar  studies 
on  other  test  sites  can  guide  the  Improvement  of  earthquake-explosion 
discrimination  and  yield  estimation. 
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